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ABSTRACT 


An experimental test of a structure results in a database of frequency response 
functions from which natural frequencies and mode shapes are identified for the given 
boundary conditions of the test. The natural frequencies of a system under a variety of 
boundary conditions can be identified by applying artificial boundary conditions (ABCs) 
at measurement locations and obtaining the frequency response function for the 
unrestrained degrees of freedom. These frequencies are found without any physical 
alterations of the test boundary conditions. Use of the ABCs in sensitivity-based model 
updating and damage detection in conjunction with baseline data provides improved error 


localization and results in a more accurate finite element model. 
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I. INTRODUCTION 


In general, a finite element (FE) model is defined by a large number of physical 
parameters, but a modal test of a structure results in a small number of modal parameters 
which are used to modify the model parameters. It is always crucial, then, to improve or 
"update" a finite element model in order to produce a reliable, confident prediction of the 
structural response of a system. Adjusting the parameters that define the model can 
reduce inaccuracies ina FE model. These parameters include, but are not limited to 
dimensional properties, moduli of elasticity, and density. Ifthe parameters in the model 
can be adjusted for closer agreement with the measured parameters (identified in a modal 
test), then the model will provide a more accurate representation of dynamic response. It 
is important to remember that no FE model can exactly predict the dynamic response ofa 
real structure, but can still provide reasonable and reliable results. This is the goal of 
model updating. 

One problem that has been historically associated with finite element model 
verification has been in measuring enough information. Procedures have been developed 
for measuring larger experimental databases. One such method that has been developed 
is known as "Perturbed Boundary Condition" testing. [Ref. 4] In this procedure, selected 
boundary conditions of a system under test are drastically changed and the system is 
retested. This method is effective, however the procedure requires that physical 
modifications be made to the structure and additional tests are necessary for each 
modification made. The result is increased time, increased work and most importantly in 


today's world, increased cost. 





A more efficient method is to identify additional and distinct mode frequencies 
from the same modal test that was performed to identify the mode frequencies of the 
standard system, without the need for physical modification of the structure. The 
additional ficaiencics and mode shapes correspond exactly to the mode frequencies 
found when combinations of measured coordinates are constrained to ground. These 
additional frequencies are therefore associated with different boundary conditions for the 
structure and no physical change in the boundary conditions has been made. [Ref. 1] 

Because the boundary conditions are not actually applied, they are termed 
"artificial boundary conditons", or "ABCs." For a single degree of freedom system the 
driving point antiresonance frequencies of a frequency response function (FRF) 
correspond to the frequencies of the structure with the driving point degree-of-freedom 
restrained to ground. The ABCs are imposed on the FE model and correspond to ideal 
constraints. With only one experimental database, this scheme yields a separate FE 
model for each ABC configuration. With the exception of the imposed boundary 
conditions, each model is identical and can be used to generate sensitivity data. In a set 
of spatially incomplete frequency response function data, the ABCs are the boundary 
conditions that define an "omitted coordinate system (OCS)," or "o-set." In performing a 
vibration test, a choice is made as to the set of coordinates to instrument with response 
transducers. This is the "analysis set,” or the "a-set." With the choice of an a-set, the o- 
set is defined as the complimentary set of coordinates. For an FE model, the o-set is of 
finite dimension, but with respect to the test system, this set is of infinite dimension. The 
natural mode frequencies of the o-set systems appear as peaks in the impedance spectra. 


[Ref. 2] 








ABC configuration frequencies are available from any set of test data due to the 
fact that a spatially incomplete FRF matrix is identically equal to the FRF matrix that is . 
calculated from the exact dynamic reduction. [Ref. 1] 

In addition to providing a larger number of frequencies for a system using one 
experimental database, the ABC also provides a means to eliminate or greatly reduce ill- 
conditioning in sensitivity equations. With the system artificially restrained at various 
measured coordinates, the columns of the sensitivity matrix which are found from the 
ABC configuration are linearly independent of the columns calculated from the baseline 
configuration. This reduces the difficulties in determining the parameters that are in 
error. Most importantly, in order to more accurately update the FE model, this 


configuration improves error localization. 








I. THEORY 


A. OMITTED-COORDINATE SET 
The equation for the steady-state harmonic response for a full order model 


(without damping) is 


k,, k,, Q? mM. m,, Xx, f, 21 

K,, k,, : mo. Mo x, 7 f, ( ) 
where Q is the frequency of harmonic excitation and f is a vector of excitations consistent 
with x and k and m are the stiffness and mass matrices, respectively. The subscript "a" 


refers to the "a-set" and the subscript "o" refers to the "o-set". Equation 2.1 is also known 


as the impedance model and can be redefined as: 

Lin Z,, x, f, 

Zea Zoo\(%o} (fo | (2.2) 
where Z(Q) is the system impedance matrix evaluated at the forcing frequency, ©. 


Assuming that {f,} = {0} (static constraint - no excitation acting on the omitted 


coordinates), the exact relationship between the o-set and a-set coordinates is as follows: 
{x,} = [1 ~ Ok:m,,| [- kk.,, z 07k, {x,} (2.3a) 


or 


{x,} =[-Z..] [Zal{x.} (2.3b) 


The origin of the o-set system results from Eq. (2.3a). By definition, the bracketed 


inverse term is 








1 


[1- 0?k='m,, jaa ~?k='m,, (2.4) 


1 
~ Det[t-0?k>'m,, 
where Det[ e ] indicates the determinant and Adj[e ] indicates the adjoint matrix. From 


Eq. (2.4), it is clear that the bracketed inverse term does not exist at those frequencies Q 
which satisfy 

Det[I- 07k3'm,,] = 0 (2.5) 
Therefore the relationship between the a-set and o-set does not exist at those frequencies. 
The frequencies that satisfy Eq. (2.5) are the eigenvalues of the system defined by [Koo] 
and [m,,], the o-set system. This system is obtained by fully constraining all coordinates 
in the a-set to ground. The use of spatially incomplete FRF data in identification of the o- 
set mode frequencies implies the use of dynamic reduction applied to the structure under 
test. [Ref. 2] 
B. EXACT DYNAMIC REDUCTION AND FRF MATRICES 

By equipping a structure with a finite number of response transducers, a reduced 
order model is defined, where the impedance of the reduced order model is not linearly 
dependent on the impedance of the full order model. The "full-order", exact FRF model 


of a structure is a FRF matrix of infinite dimension, 


H* H,, H,, 2 6 
7 H., H,, ( 


where the number of coordinates in the o-set is infinite. The FRF matrix measured in a 


test is a matrix partition extracted from the infinite dimension matrix, 


H =H (2.7) 


aa 











where the overbar notation indicates a reduced model. The matrix partition defined by 
Eq. (2.7) represents a structural dynamic model that has been reduced using exact 
dynamic reduction. From the partitioned inverse relation of the FRF matrix to its 


associated impedance, ZH = I. [Ref. 2] 

Hy, = (Ze, — Ze. ZeZon) (2.8) 
The presence of the o-set system dynamics in the term Z;) is seen in Eq. (2.8). The 
formula for the matrix inverse seen in Eq. (2.4) is applicable. Since every element in Z;) 


is singular at the natural frequencies of the o-set, Eq. (2.8) shows that elements of H 
will be singular at the o-set natural frequencies [Ref. 1]. Since the measured FRF matrix 
implicitly defines a dynamically reduced impedance model, an order n FRF matrix is 
calculated from the FE model, and the partition of this matrix which corresponds to the 
coordinates measured in a test is extracted (the a-set). The reduced model retains the 
entire modal content of the original model. 
C. DRIVING POINT FREQUENCY RESPONSE FUNCTION 

Starting again with the equation of steady state forced response for a linear 
structural dynamic system, Eq. (2.2), rewritten (with damping) as 

[k - Q’m — jOC|{x} = {f} (2.9) 
where [k -0?m - jac] = [Z] and C is the damping matrix. Transferring from physical 
to modal coordinates let 

{x} = [®]{q} (2.10) 


Rewriting Eq. (2.9), 








[Z][P]{q} = {f} (2.11) 
Premultiply by []" : 

[®]" [Z][P]{q} = [0]' {f} (2.12a) 

[Oko -27"mo + jav"Co]{q} = {3} (2.12b) 
Using orthogonality and assuming proportional damping, 

[o? - 0? + 2j6Qo]{q} = {3} (2.13) 
where «; is the natural frequency of the i mode, ¢ is the damping ratio and 
[w? ~ +2 jCQo| is a diagonal matrix. This matrix is inverted to find the modal 


frequency response: 


1 
= | ——__—__—_—__ | {5 2.14 
{4 ae | ely) 
Transforming back into physical coordinates by premultiplying by [®] and using 


{3} = []" {f}, the modal decomposition of the FRF matrix in physical coordinates is as 


follows: 


fs} 5 (9 forin (2.15) 


wo? -Q? +2j6Qa; 


where [H(Q)] = (ol i ©]'. [H(Q)] can be written: 


@? —Q? +2jC6Qa, 


ee {gr}{or}" 
[H(Q)]= > aa +2)c0e, (2.16a) 


p=l 








or any element, 


__ NRO {o°}{or} 
Hy(Q)= > o? 0" +2)C00, (2.16b) 


p=! 
This driving point FRF will be used in subsequent chapters in defining the ABC natural 


frequencies. 





10 








III. ABC CONFIGURATION FREQUENCIES 


A. DRIVING POINT ANTIRESONANCES DEFINE ABC FREQUENCIES 
FOR A SINGLE A-SET COORDINATE 
The following example using a basic 2 degree-of-freedom (DOF) system, shows 
that driving point antiresonances correspond to the natural frequencies of the structure 


with the driving point DOF constrained to ground. 


| on me 


mi ky m2 kz 


Figure 3.1 2 DOF system 


From Eq. (2.16b), the undamped driving point FRF is given by 


HQ) = See hei 


T=1 


(3.1) 


where 6; is amass normalized mode shape element, @, is the r" natural frequency, and Q 
is the forcing frequency. The frequency of the anti-resonance of H1;(Q) is given by 


le, ee ae 
oO? _ RO, +R); 


anti-res (3.2) 
j Ru+Ri 


where the modal residue is given by Rj, = 6/9;. It can be demonstrated that this 


frequency is identical to the natural frequency of the system in Figure 3.1, with the 


driving point DOF constrained to ground as in Figure 3.2. 
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ao ae = 


Figure 3.2 2 DOF system with DOF #1 restrained to ground. 


Using a simplified numerical example, if the values for Figure 3.1 are set as m)=m=1.0 


and k,=k,=1.0, the frequency of the single antiresonance 18. Qantites = Api rad/sec, which 


is identically equal to the single natural frequency of the ABC system seen in Figure 3:2; 


which is @ = 2 rad/sec. 
B. ABC CALCULATIONS FOR A 2 DOF SYSTEM 

The ABC frequencies are calculated using Eq. (2.8). A spatially incomplete FRF 
matrix is generated and inverted at each frequency. Ina plot of the elements of the 
resulting impedance matrix versus frequency, the singular frequencies are also the ABC 
system frequencies. To show an example, the 2 DOF system shown in Figure 3.1 is used 
to calculate the driving point FRF H),(Q). The system contains two modes and the 
frequencies of these modes are w = 0.618 rad/sec (f) = 0.0984 Hz) and o; = 1.618 
rad/sec (f) = 0.2575 Hz), respectively. A plot of Hi;(Q) is shown in Figure 3.3 and these 


modes can be seen on the plot. In addition to these two modes, it can also be seen that 


the anti-resonance frequency, as computed in Eq. (3.2), is Qunti-res = V2 = 1.414 rad/sec. 











0 0.2 04 O06 08 4 1.2 1.4 1.6 1.8 2 
Frequency (rad/sec) 


Figure 3.3 Hy;(Q) versus Q 
Using Eq. (2.8), the ABC system frequency can be found by calculating H;! (Q) where 
the a-set is defined as DOF #1, and H! (Q) is the scalar inverse of Hi1(Q). Hj (Q) is 
plotted versus Q in Figure 3.4. The single singular frequency corresponds to the anti- 
resonance frequency for this case (Qunti-ses = 1.414 rad/sec). This singular frequency is 
also the natural frequency of the ABC system that was obtained by constraining DOF #1 


to ground (as in Figure 3.2). Most importantly, this simple example demonstrates the 


power and usefulness of Eq. (2.8) and shows that driving point antiresonances correspond 
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3 
0 0.2 04 #06 08 1 1.2 14 1.6 1.8 2 
Frequency (rad/sec) 


Figure 3.4 Hj} (©) versus Q 
(Natural Frequency of ABC system) 


directly to the natural frequencies of the system with the driving point coordinate 
constrained to ground. [Ref. 1] 
C. ABC FREQUENCY CALCULATIONS FOR A FREE-FREE BEAM 

To further demonstrate the usefulness of the ABC technique, a free-free beam is 
examined. The model consists of 10 beam elements and each element contains two 
nodes. The nodes have 2 degrees-of-freedom, translational and rotational, giving the 
system a total of 22 DOF. The odd-numbered DOF are the translational degrees and the 
even-numbered DOF are the rotational degrees. The properties of the beam are as 


follows: length = 60 inches, EI = 5.5x10° Ibf-in’, p = 0.283 Ibf/in’, and cross-sectional 


14 





area = 0.75 in’. Response transducers (accelerometers) are placed at every other node, at 


the translational DOFs #1,5,9,13,17, and 21. This model is shown in Figure 3.5 


[\ L [\ L A L\ 


se ee 8s Ss Ss @ SG 8S  @ 8 


4 


Figure 3.5 Free-free beam - Transducer: [\ 
The set of coordinates chosen to be instrumented with response transducers is the a-set, in 
this example defined as [1 5 9 13 17 21]. It can be assumed by the given setup, that 
excitation is applied at each of the a-set DOFs, resulting in a 6 by 6 matrix. The 


impedance matrix, H 7! (Q), the scalar inverse of Hi,(Q), is calculated at excitation 


frequencies from 0-800 Hz. Figure 3.6 shows the driving point FRF for the system. 





Frequency (Hz) 


Figure 3.6 Driving Point FRF Hj;(Q) for Free-Free Beam 
(a-set DOFs: 1, 5, 9, 13, 17, 21) 
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The ABC system frequencies are found by calculating the impedance matrix, H 73 (Q). 
By plotting the 1,1 element (H 71 (Q)) of the matrix, the ABC system frequencies can be 


seen (Figure 3.7). 
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Figure 3.7 Hj) (©) for Free-Free Beam 
(a-set DOFs: 1, 5,9, 13, 17, 21) 


These ABC frequencies correspond exactly to the natural frequencies of the system with 
all of the measured coordinates (a-set) constrained to ground. This configuration is 
shown in Figure 3.8. Confirming the results, the natural frequencies of a beam with 


pin restraints in the locations shown in Figure 3.8 are calculated to be: 346.70, 385.13, 








Figure 3.8 ABC Configuration — a-set: [1 5 9 13 17 21] 
(Measured coordinates restrained to ground) 
482.73, 610.23, and 735.11 Hz. These frequencies correspond exactly to the peaks 
shown in Figure 3.7. 
A common vibrational test involves two shakers placed to excite all modes in the 
bandwidth of interest. To simulate this, the next ABC configuration set will be DOFs 1 
and 17, defined as the translations at nodes 1 and 9. Performing the same calculations as 


above, the impedance matrix, Hz} (Q), is obtained. Looking at the 1,1 element of this 
matrix, an additional set of ABC frequencies is found. A plot of Hj; (Q) over the range 


0-800 Hz is shown in Figure 3.9. The equivalent system is a beam with pin restraints at 
node 1 and node 9. The natural frequencies of this system are calculated to be 20.45, 
63.57, 112.30, 215.05, 364.61, 544.16, and 669.70 Hz. These frequencies correspond 


exactly to the peaks shown in Figure 3.9. The equivalent system is shown in Figure 3.10. 
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Figure 3.9 Hj} (Q) for Free-Free Beam 
(a-set DOFs: 1, 17) 
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Figure 3.10 ABC Configuration - a-set: [1 17] 
(Measured coordinates restrained to ground) 








IV. ABC CONFIGURATION IN SENSTIVITY-BASED UPDATING 


A. SENSITIVITY MATRIX DEFINED 

As previously mentioned, the disparity in the number of known parameters 
(measured modal parameters) versus the number of parameters to be adjusted in order to 
update a FE model defines an underdetermined problem. Another significant difficulty is 
that the parameters that are in error are unknown. Sensitivity-based updating is used for 
error localization in order to find those parameters that require adjustment. The ABC 
configuration frequencies can be used in addition to, or instead of, the standard baseline 
system mode frequencies in model updating and damage detection. The ABC 
frequencies correspond to the same structural system as the baseline frequencies, but with 
different boundary conditions. The governing equation for sensitivity-based updating is 

{Ao} = [T]{ADV} . (4.1) 
where {Ao} is a vector of natural frequency errors. The errors are the difference between 
the experimental natural frequencies and the analytical natural frequencies {@" - o*}. 
The {ADV} term is the vector of changes to be calculated for specified model 
parameters, known as "design variables," and [T] is the sensitivity matrix. This matrix is 
solved using the eigensensitivity method shown in Ref. [6]. [T] is defined as: 
Foasene} [AK] {baseline}; Where Ak = [k* — k*] and [T] is of size m x p where m is the 
number of modes and p is the number of design variables used. Each ABC system 
defines additional rows of Eq. (4.1), defining the equation, 


{Aw'} = [T']}{ADV} (4.2) 








where Ti = 00; /@DV,, and o/ is the i” natural frequency of the k" ABC configuration 


system. Baseline system quantities are identified with the superscript "0" and ABC 
systems are identified with superscripts from 1 to "k" where k is the total number of ABC 


systems used. Combining the equations, 


fao'}) {i 
oe = Mm {ADV} (4.3) 
(ao"}} (7) 


The degree of coupling between the ABC systems and the baseline system can be 
adjusted by deleting or retaining individual columns of [T*]. Partial coupling can be 
established by partitioning such that some of the DV's are associated with the baseline 
system only, some with both the baseline and the ABC system, and some with the ABC 
system only. [Ref. 1] 

Using the ABC system sensitivities helps to eliminate or greatly reduce the 
problem of poorly conditioned or rank deficient [T] matricies. Columns of [T°] can be 
replaced with columns of [T*] in order to improve conditioning. Two closely spaced 
elements in a model will always have columns of [T°] which are nearly dependent. This 
prevents the localization of the error and discrimination of the damage between the 
elements. By replacing a column of [T°] with a column of [T*], the columns are no 
longer dependent and error localization and damage detection is improved, resulting ina 


more accurately updated FE model. 
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B. IMPROVED CONDITIONING USING ABC SENSITIVITIES 

Applying the case of a free-free beam, the sensitivity matrix [T] is generated for 
the undamaged structure model in order to solve for the changes in the element EJ values 
representing damage, i.e. AEI = ADV, and 

{o°} = [T°]{AE} (4.4) 
To demonstrate, sensitivities will be calculated for all 10 elements and for the first 10 
elastic modes resulting in a [T°] matrix of size 10 by 10. A calculation of the rank of this 
matrix reveals that it is, in fact, rank deficient: Rank(T°) =5. This will not provide a 
fully determined solution for {AEI}. With excitation applied at DOFs 1 and 17, it is 
possible to identify the ABC system frequencies with these DOF constrained to ground. 
(See Figures 3.9 and 3.10) The corresponding [T’] matrix is generated and is also of size 
10 by 10, 

{o!} = [T' {AED} (4.5) 
with the calculated rank:  Rank(T’)= 10. 
The [T’] matrix is full rank due to the asymmetric ABC configuration. 
C. COMPUTER SIMULATED DAMAGE DETECTION USING ABC 
SYSTEM SENSITIVITIES 

Four examples have been constructed to demonstrate the usefulness of a single 
ABC system set of frequencies and sensitivity matrix in lieu of the standard system set. 
The ABC system described by Figures 3.9 and 3.10, a-set: [1 17], is again used. In these 
four examples, changes in elemental EI’s have been user inputted to simulate a damaged 
structure. Two subsequent examples will show the usefulness of this technique when 


damage location and magnitude is unknown. 
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1. Damage Detection: Example 1 

For this example, a 10% reduction in EI is made at elements #2 and #3. The 
natural frequencies of the undamaged (FE/Analytic) and simulated damaged 
(Test/Experimental) models are shown in Table 4.1. Using the baseline system, the 
largest full rank submatrix of [T°] is 5 by 10, and the best solution for the damage was the 


solution using modes 1 to 4. This solution is shown in Figure 4.1(a). The height of the 
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1264.6 1252.0 
1679.4 1662.9 


Table 4.1 Example 1: FE and Test Frequencies and Error 





bars in the plot indicates the magnitude of the error for each element. In this case, we 
know that the exact solution is 10% in both element #2 and element #3. If the ABC 
system sensitivities [T’] and mode frequencies are used, the largest full rank submatrix is 
10 by 10, and the best solution for the damage was the solution using all 10 modes. This 
solution is shown in Figure 4.1 (b). This was chosen as the first example because in this 
case, the baseline system provided a reliable result. It should be noted that the ABC 


system also provided reasonable results (other elements were affected by 1% or less). It 
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will be seen in the following examples that the baseline system rarely provides reliable 


results and in many circumstances actually provides results in misleading information. 





1 2 3 4 5 6 7 8 9 10 
(b) 


Figure 4.1 Damage Detection — Example 1 
(a) Baseline (b) ABC system 


2. Damage Detection: Example 2 

For this example a 10% reduction in EI is made at element #4 and a 15% 
reduction in El is made at element #6. The natural frequencies of the undamaged 
(FE/Analytic) and simulated damaged (Test) models are shown in Table 4.2. Using the 
baseline system, two solutions seemed to localize the error. The first 
was the solution using modes 1 to 5. This solution is plotted in Figure 4.2(a). The 
second was the solution using modes 1 to 4 and this solution is shown in Figure 4.2(b). 


However, it was known from the start that the damage was located at elements 4 and 6. 
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Mode FE/Analytic 
Number Hz 








Table 4.2 Example 2: FE and Test Frequencies and Error 
In a standard test, there would be no prior knowledge of damage and the results using the 
baseline system would appear to be accurate. This misleading information would not 
help in the updating of the FE model. Using the ABC system and modes 1 through 8, the 
solution shown in Figure 4.2(c) is obtained and proves to be a reliable solution based on 


the known error. 


0.2 








Figure 4.2 Damage Detection — Example 2 
(a) & (b) Baseline (c) ABC system 
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3. Damage Detection: Example 3 

For Example 3, two elements at opposite ends of the beam were chosen. A 15% 
reduction in EI was made at element #1 and a 20% reduction in EI was made at element 
#10. The natural frequencies of the undamaged and simulated damaged models are 


shown in Table 4.3. Using the baseline system, the only solution that provided any 


Mode 
Number 








10 1679.4 


Table 4.3 Example 3: FE and Test Frequencies and Error 


error localization was the solution using modes 1 to 5. This solution is shown in Figure 
4.3(a). It is immediately obvious that this solution is incorrect, as it indicates an error of 
approximately 42% in element #10 and no damage in any other elements. Using the 
ABC system and modes 1 through 7, an improved solution is obtained and is shown in 
Figure 4.3(b). 

4. Damage Detection: Example 4 

For the fourth example three elements were chosen to have reduced EI values. A 


10% reduction in E] is made at element #3, a 15% reduction at element #4, and a 20% 
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Figure 4.3 Damage Detection — Example 3 
(a) Baseline (b) ABC system 














reduction at element #5. The natural frequencies of the undamaged and simulated 
damaged models are shown in Table 4.4. The best baseline system solution was obtained 
using modes 1 to 3 and is displayed in Figure 4.4(a). The solution localized error in three 
elements, but localized damage in element #7 instead of element #4. This is obviously 
incorrect. Using the ABC system and the same modes (1 to 3), the solution displayed in 
Figure 4.4(b) was obtained. This solution also is incorrect and misleading. By using the 
ABC system and modes 1 through 6, however, an improved and correct solution is 
obtained and is shown in Figure 4.4(c). It is important to note that in this fourth example, 


the ABC system also provided one misleading solution. Given the number of modes 
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Table 4.4 Example 4: FE and Test Frequencies and Error 





used, the solution using the first six modes would be considered more reliable even if the 


damage had not been known a priori due to its improved conditioning. 
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Figure 4.4 Damage Detection — Example 4 
(a) & (b) Baseline (c) ABC system 
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For the final two examples, a random number generator is used to place errors in 
1 to 2 elements of the model by reducing EI in each of those elements by up to 50%. 
Using this method, the damage location and magnitude is unknown and left to be 
determined by using both baseline and ABC system sensitivities. Examining the element 
EI matrix from the MATLAB code checks the results determined by the sensitivities. 

5. Damage Detection: Example 5 

The damage detection results obtained by using the maximum number of modes 


possible in the baseline system , five (rank=5), are shown in Figure 4.5 (a). It appears 


1 #2 3 4 5 6 7 8 9 10 
Figure 4.5 Damage Detection — Example 5 
(a) & (b) Baseline (c) ABC system 
from this plot that element #6 is the damaged element. However, when modes 1 through 
4 are used in the baseline system (Figure 4.5(b)) it appears that #5 is the “damaged” 
element. By using modes 1 through 5 in the ABC system, Figure 4.5(c) was obtained and 
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verified that #5 was the damaged element. This was believed to be the best localization of 
error and the element EI matrix was checked for verification. The only damaged element 
was element #5 and the difference between damaged and undamaged was a 23 percent 
reduction, which corresponds to the values in Figure 4.5(b) and (c). Once again, the use 
of the ABC system provided reliable results when, without the ABC configuration, the 
baseline information alone may have been misleading. 

6. Damage Detection: Example 6 

Once again, the results using the first five modes with baseline sensitivity are 
displayed in Figure 4.6(a). Using this solution, it appears that the damage is located in 
elements #2 and #6. With the first four ads the solution shows damage in elements #2 
and #5 as shown in Figure 4.6(b). Figure 4.6(c) shows the results obtained when using 


all 10 modes and the ABC system sensitivity. It confirms that elements #2 and #5 are the 
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Figure 4.6 Damage Detection — Example 6 
(a) & (b) Baseline (c) ABC system 
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damaged elements. Checking the EI element matrix, it was confirmed that element #2 
had a 10% reduction in EI and element #5 had a 3% reduction in El. These values 
correspond to those seen in Figure 4.6(b) and (c). In summary, all 6 examples have 
provided valuable information concerning damage detection. Many times, a reasonable 
solution can be obtained using baseline sensitivities. However, with this always seems to 
be some misleading information when the baseline solutions are used. The ABC system 
will most always provide the proper solution, but is not perfect either. When the ABC 
system is used in conjunction with the baseline system, improved and accurate solutions 


are obtained. Nowhere is this clearer than in the previous two examples. 
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V. EXPERIMENTAL APPLICATION 


So far, the use of the ABC configuration to determine natural frequencies ofa 
system under a variety of boundary conditions has been demonstrated using MATLAB 
and computer generated finite elements. To fully examine this process, a laboratory 
experiment was constructed. 

A. EQUIPMENT SETUP 

A mild steel beam, 60.75 inches in length, 1.5 inches wide and 0.5 inches thick, 
was suspended from the overhead by monofilament fishing line in order to represent a 
free-free beam with two rigid body modes. The beam was instrumented with 5 PCB 
Piezotronics Inc. accelerometers, each 15.25 inches apart. Each of the accelerometers 
was wired into a 5 kHz Zonic System 7000 front end digital signal processor (DSP) 
consisting of 2 digital chassis with 16 channels each and 1 HP workstation. An excitation 
was applied successively to the opposite side of each of the accelerometers by a 25 lb 
peak force shaker and lateral excitation stand Model 2050A, manufactured by The Modal 
Shop, Inc. A MB Dynamics Model SS250VCF amplifier set to constant voltage powered 
the shaker. The shaker was operated at 0.25 volts (rms). The applied force was wired 
into the DSP as Channel 1, and the 5 accelerometers as Channels 2 through 6, 
respectively. The Zonic system was connected to the J-DEAS Master Series 5 software 


operated in the Signal Processing mode. The laboratory setup is shown in Figure 5.1. 
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B. DATA COLLECTION 

Using I-DEAS, FRFs were collected when the force was applied at each node. . 
For each reference (i.e. for the force applied at each node), the response at each node was 
measured, J-DEAS was configured to measure 1601 spectral lines over frequencies 
ranging from 0 -2000 Hz. Each shaker provided a 1601 by 5 matrix containing the 
complex amplitudes every 1.25 Hz for a given excitation. The combined result provided a 
1601 by 25 ante for the entire system. This matrix was exported to MATLAB in order 
to compute the o-set frequencies of the system. 

C. O-SET FREQUENCIES USING EXPERIMENTAL DATA 

In order to establish a comparison and validate the results obtained in the 
laboratory, a finite element model of the beam was created. The FE model contained the 
given dimensions of the beam and a Modulus of Elasticity (E) of 27.3 x 10° Ibf/in’, 
corresponding to mild steel. This E value was chosen in order to bring the baseline 
~ model frequencies closer to the test frequencies. It was composed of 64 elements with 65 
nodes (2 DOF per node for 130 DOF total). Nodes 1, 17, 33, 49, and 65 corresponded to 
transducer locations 1 through 5, respectively, on the real beam. 

The theories presented earlier and demonstrated through computer simulation 
were applied to the experimental data. By using ABC configurations, the frequencies for 
the beam under 30 different boundary conditions can be found from the baseline (free- 
free) information without having to reconfigure the experimental setup. Table 5.1 
contains the first ten natural frequencies (in Hz) of the system computed from the test and 


compares them to the natural frequencies determined by the FE model. Figure 5.2 
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contains a plot of the 1,1 element of the FRF (H(Q)) matrix and displays the first five 


modes. 
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Figure 5.2 H),(Q) for free-free beam 


34 








The first o-set system computed was the system with node | restrained to ground, 


corresponding to the system shown in Figure 5.3. The FRF matrix was inverted to find 
the impedance matrix for the a-set (Hz! (Q)) and the peaks of the impedance matrix 
correspond directly to the natural frequencies of the system with node 1 restrained to 
ground. These o-set frequencies were compared to the o-set frequencies determined by 
the FE model and are shown in Table 5.2. A plot of H;! (Q) showing the first five 


frequencies is seen in Figure 5.4. 


| a 


Figure 5.3. ABC configuration: Node | restrained to ground 
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Table 5.2 FE & Test Frequencies for Node 1 restrained to ground 
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Figure 5.4 H;!(Q) for Node 1 restrained to ground 


In both cases, the test frequencies correspond nicely with the frequencies 
predicted by the FE model. This validates the theories and expectations described in the 
earlier chapters. With only free-free data available, the natural frequencies of the system 
under any boundary configuration can be found. The number of transducers (i.e. 
measurement locations) limits the number of systems for which faces available 
and in the given example, up to 30 system frequencies can be found. Six more ABC 


systems were computed and the results for the first ten modes are tabulated in Tables 5.3 


- 5.8. The system configurations and H;'(Q) plots are shown in Figures 5.5 - 5.16. 


36 








a nll 


Figure 5.5 ABC Configuration: Nodes 1 & 2 restrained to ground 
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Table 5.3 FE & Test Frequencies for Nodes 1,2 restrained to ground 





Figure 5.6 H,! (Q) for Nodes 1,2 restrained to ground 
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Figure 5.7 ABC Configuration: Nodes 1,4 restrained to ground 
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Table 5.4 FE & Test Frequencies for Nodes 1,4 restrained to ground 
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Figure 5.8 H;!(Q) for Nodes 1,4 restrained to ground 
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Figure 5.9 ABC Configuration: Nodes 1,2,3 restrained to ground 
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Table 5.5 FE & Test Frequencies for Nodes 1,2,3 restrained to ground 
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Figure 5.10 Hz! (Q) for Nodes 1,2,3 restrained to ground 
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Figure 5.11 ABC Configuration: Nodes 1,2,3,4 restrained to ground 
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Figure 5.12 H;!(Q) for Nodes 1,2,3,4 restrained to ground 
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Figure 5.13 ABC Configuration: Nodes 1,2,3,4,5 restrained to ground 
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Table 5.7 FE & Test Frequencies for Nodes 1,2,3,4,5 restrained to ground 












[eae a 
ATT | WA TA ty 
Ni Re 
Ribas ae 
‘EA viaeiee 
PY TT te WT 
a ae 
(se ee ae Le 








Figure 5.14 H)(Q) for Nodes 1,2,3,4,5 restrained to ground 
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Figure 5.15 ABC Configuration: Nodes 1,3,5 restrained to ground 
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Figure 5.8 FE & Test Frequencies for Nodes 1,3,5 restrained to ground 
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Figure 5.16 H; (Q) for Nodes 1,2,3,4,5 restrained to ground 
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D. DAMAGE DETECTION USING EXPERIMENTAL DATA 

The next logical step to take is to try to localize error in a FE model by using the 
experimental data. In Chapter IV, sensitivity-base damage detection was used to localize 
error in a "damaged" FE model by using another FE model with no damage. We now 
have a set of data for a real beam that has no damage. The experimental data is now 
needed in order to find the error in a FE model of the same system. To perform these 
calculations, an FE model was again created, this time of 16 elements. Fewer elements 
were used in order to decrease computational time; the actual change in frequencies from 
those determined using 64 elements was negligible. A 20% reduction in EI was made at 
element #8 in the FE model. A model of the free-free beam (baseline) was created as 
well as a model of the ABC system with nodes 1 and 4 restrained to ground. 

To check the method, FE frequencies were substituted for the test frequencies in 


both the baseline system and the ABC system. The results are shown in Figure 5.17. 
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Figure 5.17 Error localization using FE frequencies 
(a) Baseline (b) ABC system 
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The error was localized using the first two modes of the baseline system (Figure 5.17(a)) 
and then was confirmed using the first two modes of the ABC system (Figure 5.17(b)). 
This validated the method and now the FE frequencies were replaced by the actual test 


frequencies for the system. These results are shown in Figure 5.18. 





Figure 5.18 Error localization using test frequencies 
(a) Baseline (b) ABC system 


The first two modes were again used for both the baseline system (a) and the 
ABC system (b) since these were the modes that resolved the error when the FE 
frequencies were used. This time, however, the baseline system localized a 20% error in 
element #8 as well as a 13% error in element #12. Clearly this is incorrect. The ABC 
system provided better results, but also indicated error in element #12 at 5%. The 


differences in the frequencies between the test data and the FE data, however slight, are 
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the sources of the discrepancies. In the first two modes of the free-free beam, the test 
frequencies differed from the FE frequencies by 0.3% and 1.0%, respectively (see Table 
5.1). For the ABC system’s first two modes, the test frequencies differed by 0.0% and 
0.4%, respectively (see Table 5.4). Therefore, in order to detect and localize errors in FE 
models with test information, it is imperative that the error between test and FE 
frequencies be minimized. Suggested approaches on how to achieve this are proposed in 


the Recommendations. 
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VI. CONCLUSIONS AND RECOMMENDATIONS 


The main objective of this thesis was to examine a method in which the natural 
frequencies of a system under a variety of boundary conditions could be found by using 
only one experimental database. Previously, experiments had to be reconfigured in order 
to find more than the baseline frequencies. This results in more time and inevitably, an 
increase in cost for structural analysis. Secondly, a method by which model error 
detection and localization could be improved was investigated. 

A. CONCLUSIONS 

The following conclusions can be drawn form the analyses presented: 

1. The natural frequencies of a structure under various boundary conditions are 
available from a baseline FRF matrix. Inverting a spatially incomplete FRF at 
each frequency results in an impedance spectrum whose peaks correspond to 
the natural frequencies of the structure restrained to ground at measure 


coordinates. 


2. With both the baseline and ABC frequencies available, structural sensitivities 
can be generated from both the baseline and from the system with boundary 
conditions artificially applied corresponding to the FRF matrix (a-set). The 
addition of the ABC data has proven to be valuable in conjunction with the 


baseline information in localizing error and/or damage. 
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B. RECOMMENDATIONS 

The following are areas of recommended improvement to continue with this field 
of study. The first two are recommendations for improving sensitivity-based error 
localization when using test data. 

1. Curve fit the FRF in -DEAS prior to exporting the files into MATLAB. 

2. Once in MATLAB, curve fit the impedance spectra near the ABC frequency 

to identify the ABC natural frequencies. 
3. Continue with the use of ABC systems on more advanced structures to further 


this study. 
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APPENDIX MATLAB CODE 


A. MAIN PROGRAMS 

Build2Beams.m 

% This program assembles the mass and stiffness matrices for 2 beams, with user 

% defined boundary conditions referred to as "BeamA" (analysis) and "BeamX" 

% (experimental). The program can be run in both “build” and “read” modes. In the 
% “build mode,” the user provides baseline data for BeamA, assumed to be a 

% homogeneous, uniform beam. Data provided: 

% 

% (1) Beam length 

% (2) Number of elements 

% (3) Nominal EI 

% (4) Nominal cross-sectional area 

% (5) Nominal weight density 

% 

% The program then prompts the user for instructions on how to modify "BeamA" 
% data to arrive at "BeamX" data. The user can modify element EI values and the 
% modification can be applied to either a single element, or range of elements. The 
% beam definition data is saved in a binary (.mat) file "beamdata" at the end of 

% execution. 

% 

% In the "read" mode, the program will load the beam data written at the last 

% execution of the program and will then allow the user to modify the data for 

% BeamX. By answering "y" to this prompt, the program will copy BeamA data into 
% BeamX data, and then prompt the user (identically as described above) for data on 
% how to adjust the BeamA data in order to arrive at the BeamX data. 

% 

% The program calls the following support programs: 

% 

% BeamA_Prompt.m -- Prompts User for BeamA nominal beam data 

% (“build” mode only) 

% BeamX_Prompt.m -- Prompts User for BeamX modification beam data 

% BoundaryConditions.m -- Prompts user for B.C.'s and applies them. 

% Assemble2Beams.m _ -- Called by Build2Beams, builds [ka] [ma] [kx] [mx], 

% plots freqs. 

% Saves data to "beamdata.mat" 

% clearspace.m -- Clears all unnecessary variables from workspace. 

% 

% Written by J.-H. Gordis 

% Modified by C.P. DeGregory 
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disp(''); 
disp( ' Build2Beams execution...'); 


input_flag = 'b’; 
input_flag = input(' "B"uild new beams or "R"ead from file ? (b/r) ','s}); 


% Build_or_read if loop: 


OY ESE C ORC EFE EE SIAEISIOIOIOK SSSI ICICI IAC IEA AICI ICI RCAC ICA FOR ICR AAAI A tek 


if input_flag =='b'; % (IF#1) Prompt user for beam definition data 
Of EEE SS ICIS IEE ISIS CEE REIS ORICIC ICICI AACA AICI AA Aa A AA 


disp(' Building 2 beams from scratch...’) 


BeamA Prompt;  % Prompt for BeamA Data: run prompt script 

BeamX Prompt;  % Prompt for BeamX Modification Data: run prompt script 
Assemble2Beams; % Run script to assemble mass and stiffness. 
BoundaryConditions; % Prompt for, and apply boundary conditions 


% *** Save Defining Parameters for Beams *** 


disp(' ...saving beam data to file’) 

save beamdata connect ndof element_length element_EI element_mass mass_diag 
restraint_switch free_dof_set... 

num_elements num_rbm lama lamx 

%save beamdata phia phix 


OER A TCI IG A AC I TOI ICI I ICICI ka a a a a ae a a a ae sk a CAR A aC ae IA A He A 


elseif input_flag =='r'; % (IF#1) file read for beam definition data 


OIE AACA CHIE ISSAC I RGIS AACR ACCIACCA A A ICR ACI ICH TCI AA A A A 


disp(' Loading 2 beams from ascii file "beamdata"...’) 
load beamdata 


% Prompt user to either build two beams from data in file, or modify beamx data: 
modify beamx ='n’; 
modify _beamx = input(' Modify beamX data (y/n)? ','s’); 


if modify _beamx =='y'; % (IF#2) 
disp(’ Modify BeamX...’) 
reset_beamx ='y'; 
reset_beamx = input(' Reset BeamX properties to BeamA properties? ... 
(y/n) ','s'); 
if reset_beamx =='y'; % Copy beamA data to beamX data: 
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element_FI(:,2)=element_EI(:,1); 
element_mass(:,2)=element_mass(:,1); 
% ...Display BeamA properties to user: 
sprintf((Number of elements: %3i',num_elements) 
end; % end reset_beamx 


BeamX_Prompt; % Prompt for BeamX Modification Data: run script 


% End modify_beam if block (IF#2) 
end 


Assemble2Beams; % Run script to assemble mass and stiffness. 
BoundaryConditions; % Apply boundary conditions (that already exist) 


% *** Save Defining Parameters for Beams *** 


disp(' ...saving beam data to file’) 
save beamdata num_elements connect ndof element_length element_E] ... 
element_mass mass_diag restraint_switch free_dof_set... 
lama phia lamx phix num_elements connect ka ma kx mx num_rbm 


Of) FACES SS IG ICRI IE ACI IAC A ICI AA CI I ACR I A CE ACR a ACR ICC Eo 


end; % End Build or read if loop IF#1 


OES REGS ER IGG ERAGE RCC IOIC I ACCA II ACI AGE SF 0% a 05 A A a A 


disp(' Build2Beams end.') 
% 
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Beam_Sensitivity.m 


% This program calculates the stiffness and mass matrix partials by finite 

% differences. That is, for example, the [k] matrix is assembled twice, once for the 
% nominal beam parameters, and a second matrix is assembled based on a small 

% perturbation of element EI. 

% 


% This program makes use of the beam data created by the program 
% “Build2Beams.m,” reads the data file BeamData and resets BeamX EI data to be 
% the same as BeamA data. 


% 

% The user is prompted for element label numbers to be used as design variables by 
% the script “BeamSens_ Prompt” which stores this element label data in the matrix 
% “EI el_lbl.” 

% 

% El_el_Ibl: This array is (max_num_dv,num_elements), where max_num_dv is 
% the maximum number of design variables allowed, and each design variable may 
% include from 1 to "num_elements" elements. Therefore, the sets of elements 

% included in two design variables may share common elements. 

% 

% Based on the number of non-zero rows in the array “EI_el_lbl,” matrix partials 
% are calculated for each non-zero row. For example, to start, the first row of 

% EI el_Ibl is read, and a small perturbation in element EI is applied simultaneously 
% to the elements listed in this row. 

% 

% The script generates a matrix containing mode frequency sensitivities: 

% 

% [sens EI]: num_modes * size(EI_el_Ibl,1) 

% 

% 

% Written by J.H. Gordis 

% Modified by C.P. DeGregory 
% 


OSI IIR ICI ICICI ICI ICICI ICCC AGI ARR I ECC A A AR 


% Start code: 
ae 





format long; 


disp("'), 
disp(’ BeamSensitivity execution...');disp(’ '); 
disp(' Loading 2 beams from ascii file "beamdata"...');disp(’ '); 
load beamdata 
num_elements = num_elements 
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clear num_mass_ dv num_EI_dv define M_K_ dv EI_el_lbl mass_el_Ibl 
clear sens EI sens_mass 


BeamSens Prompt; % Prompt for DV element label definitions 


% 24 2G 2 she fe 2K Ke fe ah 2 2 oe fe fe a 2 oft 2 fe ofc eis fe 2 CC oe of fe a oe 2 eo 2 2 oe of ee oe 2 2 oo oi fe ke oe 2 ee 2 2 2k 2k 2 2 ok ok 
% 2 oe 2 2 2c 2s fs fs os ofc ae oft fe fe ale ale afe ofe ofc oie 2s fe fe afe oe INITIALIZATION oe fe ok 2 oe eo oe fe ok ke oe ok fe ok fe ae 2 fe 2k 2k oe 2 2k 2k ok oe 
% Ke fe fe ke 2 fe ak 2 2k oi fe ae oe ke oe 2 ie oe 2 oe ke oe 2c ee 2 oe ie fe 2 2 2 oe fe 2A 2 fe oe 2c fe 242 eg 2s 2 ke 2g 2 ae 2 2 2 2 0k 2 2 oe 2K Ok 


mass change=1; % Percentage mass change for sensitivity calculation. 
EI change = 1; % Percentage EI change for sensitivity calculation. 
element_mass_orig = element_mass; % Copy properties over to retain them. 
element_EJ orig = element_EI; % 


icnt_ dv =0; % Initialize number of DV's to zero. Will count up all DV's 


% Prompt for number of mode frequencies to generate sensitivities for: 
fm 





[LL tt OL Lt Dt dd DOLLS OS PDD OSD 





LLL lt Od 





ene 





num_modes = input(' Enter number of elastic modes for sensitivity calculations>> '); 
start mode = num _rbm+1; % Skip the rigid body modes. 


% fe fe 2h fe of fe fe ok of 2 fe fe ake ae ke of nc 2c oe of oe fe fe 2 fe 2 fe se ake fe ake ae ke ae fe fe oe fe se fe ok oie fe fe ake oe ok 2 ke ake fe 2k oie fe ois 9k ke oie oe oie 2 ie 2k 2k ok ok 2K ok 


% 2fe fe fe fe fe ke of fe of ac ofc 2k 2k fe oe oe oe fe oe of ake ofc 2 ofc it 2s a 2 fee af fe fe oe fe ae oe ac oe fe ie 2 fe 2K fe 2k 2 oe 2 2 he ae of oe ie 2s 2 oe oi 2k oe ok 2 OK 
num_EI_ dv = size(EL_el_lbl,1); 
for icnt_dv = 1:num_EI_dv; % Repeat all of this for all EI DV's defined. 


% Resetting BeamX properties to BeamA properties...');disp(’ '); 
element_EI(:,2) = element_EI(:,1); 


% Get element labels for current DV: 
Ibls = EI_el_Ibl(icnt_dv,find(EI_el_Ibl(icnt_dv,:))); 


% Perturb EI props for all elements in "bls": 
element_EI(Ibls,2) = element_EI(Ibls,2) * (1 + (EI_change/100) ); 


Assemble2Beams; % Run script to assemble baseline and perturbed beams. 
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end; 





BoundaryConditions; % Apply boundary conditions. 


[lam_base,phi_base]=fmodes2(ka,ma); % Get modes of baseline beam 
% Form EI derivative matrices: 
k_delta = (kx - ka)(EI_change/100); % in %/100 


% Mode freq sens loop: 
end_mode = start_mode + (num_modes - 1); % Get requested number of modes. 
row_num = 0; 
for icnt_modes = start_mode:end_mode; 
row_num = row_num +1; 
sens El(row_num,icnt_dv) = phi_base(:,icnt_modes)' *... 
k_delta * phi_base(:,icnt_modes); 
end; 


% End "for icnt_dv" outer loop for sensitivity calculations 
% Run script to clear out variables: 


clearspace 


% Copy element EI and mass properties back into arrays: 
element_EI =element_EI_orig; 
element_mass = element_mass_orig; 


save beamdata element_EI element_mass sens_EI 
% Save all variables to beamdata (essentially does an append to existing beamdata) 


disp(’ BeamSensitivity end.') 


% 
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Beamex.m 


This program calculates an FRF matrix and the inverse of an FRF matrix in order 
to identify the Artificial Boundary Condition frequencies. It loads the binary 
(.mat) file “Beamdata” and determines the natural frequencies of the system from 
the stiffness and mass matricies by calling the function “fmodes.” The program 
also calls the function “fOset_from_Aset” to compute the o-set from a user 
defined a-set, as well as “fmodes” again, this time to solve for the o-set natural 
frequencies (which correspond to the ABC frequencies). The FRF is 

calculated as the inverse of the Impedance (H = Z-1, where Z =k - (2m). Both 
the FRF and inverse FRF (Impedance) plots are displayed over the frequency 
range 0 — 800 Hz. 


load beamdata 


[lam,phi] = fmodes(kx,mx); 
freq = sqrt(lam)/2/p1; 


% 


User is asked to input a-set 


aset = input('Enter a-set (vector form): '); 


[oset] = fOset_from_Aset(length(mx),aset); 


ko = 


kx(oset,oset); 


_mo = mx(oset, oset); 


disp('‘A-Set Defined as:') 


aset 


disp('O-Set Frequencies’) 
[lamo,phio] = fmodes(ko,mo); 
freqo = sqrt(lamo)/2/pi; 

icnt = 0; 

freq = [1 : 1 : 5030}; 

freqp = ifreq / 2 /pi; 


for ifreq = ifreq; 


sia =icent+ 1; 
= inv(kx + 001*sqrt(- 1) * kx- ifreq"2 * mx); 
ae h(aset,aset); 
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hxp1 1(icnt) = hx(1,1); 

% hxp22(icnt) = hx(2,2); 
% hxp33(icnt) = hx(3,3); 
% hxp44(icnt) = hx(4,4); 
% hxp55(icnt) = hx(5,5); 
% hxp66(icnt) = hx(6,6); 


zabc = inv(hx); 

zabc11(icnt) = zabe(1,1); 

% zabc22(icnt) = zabc(2,2); 
% zabc33(icnt) = zabc(3,3); 
% zabc44(icnt) = zabc(4,4); 
% zabe55(icnt) = zabce(5,5); 
% zabce66(icnt) = zabc(6,6); 


end 


% Plot FRF and Impedance 
clf 

figure(1) 
plot(freqp,log10(hxp11));grid 
xlabel('Frequency (Hz)’) 
axis([0,800,-7,3]) 


figure(gcf+1) 
plot(freqp,log10(zabe11));grid 
xlabel('Frequency (Hz)') 
axis([0,800,-1,7]) 


save hdata freqp zabe1 1 
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Solveos.m 


% This program solves for the vector {(EI} following user guidelines. It first loads 
% data which was saved after two runs of the “Build2Beams” and 

% “Beam_Sensitivity” programs. The first run is fora baseline free-free beam with 
% given element EI modifications. The analytical and experimental (simulated 

% damaged) eigenvalues (lama & lamx) and the sensitivity matrix (sens_ET) are 

% stored in beam0.mat. The second run is for a beam with user defined constraints 
% at given measurement locations (a-set) and lama, lamx and sens_EI for this case 
% are stored in beam1.mat. 

% 

% Written by J.H. Gordis 

% Modified by C.P. DeGregory 
load beam0 

lama0=lama; 

lamx0=lamx; 

TO =sens_EI; 

clear lama lamx sens_EI 

load beam 1 

lamal=lama; 

lamx1=lamx; 

T1 =sens_E]; 


clear lama lamx sens_EI 


rerun ='y'; 
while rerun == 'y' 


% Prompt for number of Beam mode frequencies: 


Os ei a 


Note: for free-free structures, only elastic modes are in T matrix 


% 





[LLL D at td 





num_modes0 = size(T0,1); 
num_modesi = size(T1,1); 


sprintf('‘Total Beam0 mode sens. data avail: %3i', num_modes0) 
sprintf(‘Total Beam1 mode sens. data avail: %3i', num_modes1) 


disp(''); disp(’ For Beam0: "); 
use_all_ modes = input(' Use all Beam0 modes ? (y/n) >> ','s’); 


if use_all_ modes == 'n'; 
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disp(' Enter (elastic) mode numbers to include...) 
disp(' Use MATLAB vector format> 13 5:79 ') 
use_modes0 = input(' >> ','s’); 
use_modes0 = eval(['[',use_modes0,']']); 

else 
use_modes0 = 1:num_modes0; 

end: 


disp(''); disp( For Beam1: '); 
use_all_ modes = input(' Use all Beam! modes ? (y/n) >> ','s'); 


if use_all_modes == 'n’; 
disp(' Enter (elastic) mode numbers to include...") 
disp(' Use MATLAB vector format> 13 5:79 ') 
use_modes1 = input(' >> ','s’); 
use_modes] = eval({'[',use_modes1,']']); 

else 
use_modes! = 1:num_modes1; 

end; 


% Prompt for number of DV's: 


rm 








eee noeat eee aeNer Say 


num_DV0 = size(T0,2); 
num_DV1 = size(T1,2); 


sprintf(‘Total Beam0 DVs avail: %3i', num_DV0) 
sprintf(‘Total Beam! DVs avail: %3i', num_DV1) 


% Prompt user for number of DVs to include: 


disp(''); disp(' For Beam0: '); 


sprintf(‘Number of Beam0 EI DVs: %3i',num_DV0) 
use_all_dv = input(' Use all Beam0 DVs ? (y/n) >> ','s’); 
ifuse_all dv ~'y’; 
disp(' Enter DV label(s) to include...') 
disp(' Use MATLAB vector format> 13 5:79 ') 
use_dv0 = input(' >> ','s'); 
, _use_dv0 = eval(['['juse_dv0,']']); 
else 
use_dv0 = 1:num_DV0; 

end 
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disp(''); disp(' For Beam1: '); 
sprintf(‘Number of Beam1 EI DVs: %3i',num_DV1) 
use_all dv = input(' Use all Beam1 DVs ? (y/n) >> ','s'); 
if use_all_ dv ~='y’; 
disp(' Enter DV label(s) to include...") 
disp(' Use MATLAB vector format> 13 5:79 ") 
use_dvl =input(' >>','s’); 
use_dv1 = eval(['[',use_dv1,']']); 
else 
use_dvl = l:num_DV1; 
end : 


% Assemble Composite Sensitivity Matrix: 


es ERE tet set RP oD tee ae Ee ea eee, 


format short 


Tboth = [TO(use_modes0,use_dv0);T1(use_modes1,use_dv1)]; 


% Calculate Frequency Error Vector for the modes retained: 


if a pe Seta SIN reat it aS 





(Ot tt tt 


del_freq = [lamx0(use_modes0)-lama0(use_modes0);lamx1(use_modes1)-... 
lamal(use_modes1)]; 


% Solve system for DV changes: 
% WX 








a 





LPL LN OP Id Ot OH me 


sprintf("Min Dim of T: %3i',min(size(Tboth))) 
sprintf(‘Rank of T: %31',rank(Tboth)) 
sprintf('‘Condition of T: %5.3e',cond(Tboth)) 
delta_dv = Thoth \ del_freq 
bar(1:10,abs(delta_dv)),grid 


rerun = input(' Rerun for different case ? (y/n) >> ','s'); 


end; % end rerun while loop 
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Main.m 


% This program is similar to the “beamex” program in that it computes FRF and 

% Impedance matrices. This script, however, is for a spring-mass system. In 

% particular, this program solves for the system shown in Figure 3.1. The stiffness 
% and mass matrices are computed using the function “fspringmass” and the values 
% for connectivity, mass specification and fixity. The values required for the system 
% shown, in addition to the a-set (DOF #1), are contained in the script “testdata.” 

% The 1,1 elements of the FRF and Impedance matrices are plotted versus 

% frequency. 

clear 

testdata 


[k,m] = fspringmass(conn,mspec,fix); 
[lam,phi] = fmodes(k,m); 


[oset] = fOset_from_Aset(length(m),aset); 


ko = k(oset,oset); 
mo = m(oset,oset); 


icnt = 0; 
ifreq = [0:.01:4]; 
freqp = ifreq/2/pi; 
h = zeros(size(k)); 
for ifreq = ifreq 
icnt =icnt + 1; 
h = inv(k — ifreq*2*m); 


hxp1 1(icnt) = hc1,1); 

%hxp22(ient) = h(2,2); 
%hxp33(icnt) = h(3,3); 
%hxp44(icnt) = h(4,4); 
%hxp55(icnt) = h(5,5); 
%hxp66(icnt) = h(6,6); 


zabe11(icnt) = inv(h(1,1)); 

%zabc22(icnt) = inv(h(2,2)); 
%zabce33(icnt) = inv(h(3,3)); 
%zabc44(icnt) = inv(h(4,4)); 
%zabce55(icnt) = inv(h(5,5)); 
%zabc66(icnt) = inv(h(6,6)); 

end 
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orient landscape 

figure(1) 
plot(freq*2*pi,log]10(hxp11)), grid 
xlabel(‘Frequency (rad/sec)’) 
axis([0,2,-2,3]) 


figure(gcf+1) 
plot(freq*2*pi,log10(zabe11)), grid 
xlabel(‘Frequency(rad/sec)’) 
axis([0,2,-3,2]) 
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Hmat.m 


% This program is used for the experimental data obtained using the I-DEAS 


% software. It takes the data stored in beamb and creates an FRF matrix as well as 
% an inverse FRF matrix for any boundary condition. The user inputs the 

% frequency range and the nodes to be used in the a-set. 

% 

% Written by C.P. DeGregory 

clear 


load beamb % Loads data file containing complex amplitudes for each frequency data 
% point in the spectrum. 

hl1 = Hbase(:,1); 

% Plot the FRF for the free-free beam (There are 1601 spectral lines: one every 1.25 Hz) 


plot(0:1.25:2000,20*loag1 0(abs(h1 1))),grid; 
axis([{0,500,-40,50]); 


first = input(‘Enter first frequency: ‘); 
last = input(‘Enter last frequency: *); 


start = (first/1.25) + 1; 
finish = (last/1.25) +1; 


% Initialize 

icnt = 0; 

H = zeros(5*((finish-start)+1),5); 

freqp = [first:1.25:last]; 

Zap = zeros(length(freqp), 1); 

aset = input(‘Enter nodes to be used in a-set (vector form): ‘); 
% Build H matrix 


for i = (start) : (finish) 


H((5*ient-4):5*icnt,:) = [Hbase(i,1:5); Hbase(i,6:10); Hbase(i,11:15);... 
Hbase(i, 16:20); Hbase(i,21:25)]; 


% Build Ha, matrix 


Ha((length(aset)*icnt)-length(aset)+1: (length(aset)*icnt), : ) = H ((5*... 
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icnt-1) +aset, aset); 

Za((length(aset)*icnt)-length(aset)+1: (length(aset)*icnt), : ) = inv(Ha((length... 
(aset)*icnt)-length(aset)+1: (length(aset)*icnt), : ); 

% 1,1 Element for plotting 

Zap(icnt) = Za(((length(aset)*icnt) — length(aset)) +1,1); 


end 


plot(freqp, 20*log10(abs(Zap))),grid 
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B. SUPPORT PROGRAMS 


Assemble2Beams.m 


% This program assembles the global mass and stiffness matrices of the two beams 
% defined in “Build2Beams” by combining element k and m matrices solved in the 
% function “fbeamkm.” 

% 

% Loop over the two beams: 


for icnt_beams = 1:2; 


k=zeros(ndof,ndof); 
m=zeros(ndof,ndof); 


% Loop over the number of elements in the structure: 


PL LR Lt I ad dd 





OF eet assets Se tad ae ae ca aie 


for elnum=1:num_elements; 


dofl=2*connect(elnum, 1)-1; 
dof2=dofl+1; 
dof3=2*connect(elnum,2)-1; 
dof4=dof3+1; 


% ... set up beamel function call: 


1 temp=element_length; % Using fixed element lengths 
ei_temp=element_EI(elnum,icnt_beams); 
m_temp=element_mass(elnum,icnt_beams); 


[kbeam,mbeam]=fbeamkm(l_temp,ei_temp,m_temp); 


k(dofl :dof2,dofl :dof2)=k(dofl :dof2,dofl :dof2)+kbeam(1:2,1:2); 
k(dofl :dof2,dof3 :dof4)=k(dofl :dof2,dof3 :dof4)+kbeam(1:2,3:4); 
k(dof3 :dof4,dofl :dof2)=k(dof3:dof4,dofl :dof2)+kbeam(3:4, 1:2); 
k(dof3 :dof4,dof3 :dof4)=k(dof3:dof4 dof3:dof4)+kbeam(3:4,3:4); 


m(dof] :dof2,dofl :dof2)=m(dofl :dof2,dofl :dof2)+mbeam(1:2,1:2); 
m(dofl :dof2,dof3 :dof4)=m (dof! :dof2,dof3:dof4)+mbeam(1 :2,3:4); 
m(dof3:dof4,dofl :dof2)=m(dof3 :dof4,dofl :dof2)+mbeam(3:4, 1:2); 
m(dof3:dof4 ,dof3:dof4)=m(dof3 :dof4,dof3 :dof4)+mbeam(3:4,3:4); 


% end loop over the number of elements: 
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end 


% Reassign k and m to new variables and add lumped masses 


yarn Spee ee Srey eee eRe eee Nese 6 ee re 


if icnt_beams == 1; 
ka = k;ma = m; 
elseif icnt_beams == 2; 
kx = k:mx = m; 

end 


% End icnt_beams loop: 
end 


clear dofl dof2 dof3 dof4 1_temp ei_temp m_temp icnt_beams 
clear k m kbeam mbeam elnum 
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BeamA_Prompt.m & BeamProperties.m 


% 
% 
% 
% 


This program prompts the user for Beam A data. The first inputs required are the 
nominal physical properties. The short program file “BeamProperties” is offered 
to the user. The data listed in “BeamProperties” is the data that was used for all 
beams throughout this thesis. 


disp("');disp(’ '); 
disp(' Enter nominal physical properties for first beam’) 


displ eee 








t 
a ea ie 


props file ='n'; 
props_file = input(' Run "BeamProperties.m" script to define nominal properties ?... 


(yin) >> ','s'); 


if props_file =='y’; 


else; 


end; 


BeamProperties; 


total length =input(' Enter the total beam length in inches: '); 

num elements = input(' Enter the number of elements: '); 

nominal EI = input(' Enter the nominal EI value (Ibf/in“2): '); 
nominal area =input(' Enter the nominal cross section area (in%2): '); 
nominal density = input(' Enter the nominal weight density (Ibf/in%3): '); 


for icnt_elements = 1:num_elements;, 


end 


connect(icnt_elements,1:2) = [icnt_elements,icnt_elements+1]; 


ndof=length(connect)*2+2; % this is unrestrained beam DOF 


element_length = total_length/num_elements, 

element_EI(:,1) = nominal_EI * ones(num_elements,1); 

element_area(:,1) = nominal_area * ones(num_elements,1); 
element_density(:,1) = nominal_density * ones(num_elements,1); 
element_mass(:,1) = (element_density .* element_area * element_length)/386.4; 


% End Data input for 1st beam: Copy data for 2nd beam: 


Op ta eae Nore sa ao, 








RE I 


element_EI(:,2) = element_EI(:,1); 
element_area(:,2) = element_area(:,1); 
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element_density(:,2) = element_density(:,1); 
element_mass(:,2) = element_mass(:,1); 


% End BeamA_Prompt.m 


% 





% 





% BeamProperties.m 
% This is the file to load nominal beam data. 


total length = 60; 
num_elements = 10; 
nominal EJ =5.5e5; 
nominal area =.75; 
nominal_density = .283; 
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BeamX_Prompt.m 


% This program prompts the user for Beam X data. Specifically, the user is asked to 


% choose whether to input element EI modifications or to have them chosen 

% randomly. If chosen randomly, the data is saved in the binary (.mat) file 

% “changedata” and is then offered as a choice for the next run of this code. By 

% doing this, the option exists to create 2 free-free beams with unknown damage and 
% 2 a-set constrained beams with the same (but still unknown) damage. 


disp(' ');disp('’); 
disp(' Modify nominal physical properties for second beam’) 


disp(’ ‘Laie ides sped pect a psec 8 paps aCe ie ets tsi en DR at IAL preg rosea ~~~) 


% Adjust EI values for second beam: 


option = 'u'; 
option] ='n'; 
change EI ='n’; 


change EI = input(' Modify single/range element EI values (y/n)? ','s’); 
option = input(’ User defined, random, or neither (u,r,n)? ','s'); 
option] = input(' Read previous modification values (y/n)? ',’s'); 
if option! =='y’ 
load changedata 
change EI ='n' 
end 


while (change EI ~='n') & (option == 'u’) 
disp(' Enter element label(s) for EI modification’) 
disp(’ Use MATLAB vector format> 1 3 5:79 ') 


new_lbls = input(’ >>','s'); 
new_lIbls = eval(['['snew_lbls,']']); % Converts string to vector SE labels 


disp(' Enter EJ change for element range’) 

EI_change = input(' Enter percentage E] change (+/- %) '); 
element_EI(new_Ibls,2) = element_EI(new_Ibls,2)+... 
(EIl_change/100) * element_EI(new_Ibls,2); 

change EI = input(' Modify another element EI value (y/n)? ','s’); 


end; % end while 
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if option == 'r' 

numels = fix(3*rand+1); 

for i=1:numels 
new_lbls = fix(10*rand+1); 
format short 
E]_ change = -rand/2; 
element_EI(new_Ibls,2) = element_EI(new_Ibls,2)+-.. 
EI_change*element_EI(new_Ibls,2); 

end 


save changedata numels new_lIbls element_EI 


end 


% End BeamX_Prompt.m 
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BeamSens_Prompt.m 


% This program prompts the user for design variable element label defintitions. It is 
% called by the program “Beam_Sensitivity.” 


max _num_ dv =50; 

% ote of oi ok oie ke 2h ok ok oie ok 2k 2k ok ok EI DV ELEMENT LABEL PROMPTING oh ok 2 fe oe oe ke ok 2 ok ok oe ok oe ok 
define EI_dv =input(' Calculate element EI sensitivities (y/n)? ','s'); 

if define EI dv =—'y'; 


all_elements = input(' Define EI DV"s for all elements? (y/n) ’,'s'); 
if all_ elements == 'y’; 
for icnt_elements = 1:num_elements; 
EI_el_lbl(icnt_elements,1) = icnt_elements; 
end; 
define_E]_ dv ='n’; 


else; % all elements == no, user selects elements 
for icnt_dv = 1 : max_num_dv; 


sprintf(' Enter element label(s) for EI DV # %3i',icnt_dv) 
disp(' Use MATLAB vector format> 1 3 5:7 9 -- "d" to finish’) 
new_lIbls = input(’ >> ','s’); 
if new_lbls ~= 'd'; 
new_Ibls = eval(['[,new_Ibls,']']); % Converts string to 
% vector of labels 
EI _el_Ibl(icnt_dv,1:length(new_Ibis)) = new_Ibls; 
else 
break 
end; 


end; % for icnt_dv 
end; % if all_elements 


mass el lbl=[ ]; 
end; 


% End script. 
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Boundary_Conditions.m 


% 
% 
% 
% 
% 


This program prompts the user for boundary condition information, creates a 
vector of DOF (with respect to the unrestrained structure) and then extracts the 
rows and columns of the complementary DOF. The program defines vector 
"free_dof_set" containing the list of unrestrained DOF. The boundary conditions 
for the given system are applied in this script. 


if exist(‘free: dof_set'}==0; % Build free_dof_set vector 


disp(' Select a boundary condition set:') 
disp(' (1) Clamped-free’) 

disp(' (2) Clamped-Clamped') 

disp(' (3) Pinned-Pinned’) 

disp(' (4) User-Defined’) 

disp(' (5) Free-Free') 


BC_Choice = input(' >> Enter choice: '); 
if BC_Choice==1; % Clamped-free 


free_dof_set = [3:ndof]; 
restraint_switch ='y'; 


elseif BC_Choice == 2; % Clamped-Clamped 


free_dof_set = [3:ndof-2]; 
restraint_switch = 'y’; 


elseif BC_Choice == 3; % Pinned-Pinned 


free_dof_set = [2:ndof-2 ndof)]; 
restraint_switch = 'y'; 


elseif BC_Choice == 4; % User-Defined 


icnt_dof = 0; 
add_dof='y'; 
while add_dof =='y'; 
be_node = input(' Node number for restraint ? "0" to end: '); 
if bc_node == 0; 
break 
end; 
be_coord = input(' Translation or Rotation ? (t/r) ','s'); 
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icnt_dof=icnt_dof+ 1; 
if bc_coord == 't’; 

be_DOF(icnt_dof) = 2 * be_node - 1; 
elseif mass_coord == 'r'; 

be_DOF(icnt_dof) = 2 * be_node; 
end; % End if-elseif block 


end; % End while add_dof 
bc_boolean = ones(ndof,1); %[111... icnt_dof] 
be_boolean(be_DOF) = zeros(length(bc_DOF),1); % Put zeros in restrained dof 
all_dofs = [1:ndof]; % List of all dof 


free_dof_set = all_dofs.*bc_boolean’; % Extract free dof 
restraint_switch = 'y’; 


elseif BC_Choice == 5; % Free-free beam 


free_dof_set = [1:ndof]; 
restraint_switch = 'n’; 


end; % End if-elseif choice block 
end; % End exist block 
free_dof_set(free_dof_set==0)=[ ]; 
ka = ka(free_dof_set,free_dof_set); 
ma = ma(free_dof_set,free_dof_set); 
kx = kx(free_dof_set,free_dof_set); 


mx = mx(free_dof_set,free_dof_set); 


% 3 3 2c fe fe 2 ie ck 26 ok ok oe ok Ok END BOUNDARY CONDITIONS.M 2 2 2 eo oe 2k 2c oie ake ok oe 
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C. FUNCTIONS CALLED BY MAIN & SUPPORT PROGRAMS 


function[{kbeam,mbeam] = fbeamkm(1_temp,ei_temp,m_temp) 


% This function computes the elemental stiffness and mass matrices defined using 
% Galerkin’s Method and Hermitian shape functions. [Ref 7] 


krowl = [12 6*1_temp -12 6*1_temp]; 

krow2 = [6*1_temp 4*]_temp”2 -6*]_temp 2*1_temp%2]; 
krow3 = [-12 -6*1 temp 12 -6*|_temp]; 

krow4 = [6*1_ temp 2*]_temp”2 -6*]_temp 4*1_temp*2]; 
kbeam = (ei_temp/l_temp%3)*[krow];krow2;krow3;krow4]; 


mrow] = [156 22*1_ temp 54 -13*l_temp]; 

mrow2 = [22*1_ temp 4*1_temp’2 13*1_temp -3*1_temp”2]; 
mrow3 = [54 13*1 temp 156 -22*1_temp]; 

mrow4 = [-13*1_temp -3*1_temp*2 -22*]_ temp 4*1_temp”2]; 
mbeam = (m_temp/420)*[mrow1;mrow2;mrow3;mrow4]; 


function [phi_normal,orth]=massnormal(phi,mass) 
% This function mass normalizes a modal matrix. 


a=size(phi); 

nummodes=a(1,2); 

phi_normal=zeros(size(phi)); 

% 

- for i=] snummodes; 
modalmass(ii)=phi(:,ii)'*mass* phi(:,11); 
if modalmass(ii)~=0, 

phi_normal(:,ii)=(1/sqrt(modalmass(ii)))*phi(:,11); 
else 

phi_normal(:,ii)=phi(:,ii); 
end 

end 

% 

% do the orthogonality calculation: 

% 

orth=phi_normal'*mass*phi_normal* 100; 
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function [omega,phi]=fmodes(k,m) 


% This function returns a vector containing eigenvalues (rad/ sec)’ and a matrix 
% containing the mass normalized mode shapes. The mode information is sorted by 
% frequency in ascending order. 


[v.d]=eig(m\k), 
[omega,index]=sort(diag(d)); 


[phi,orth]=fmassnormal(v(:,index),m); 


format long 

disp(‘Freqs in Hz.:’) 

disp( ~~~ ') 

disp(sqrt(omega(1 :min(length(omega),12)))/2/pi) 
disp('------------- ‘ 


function [oset] = fOset_from_Aset(ndof,aset); 


% This function determines the complementary o-set from a set [1:1:ndof] and the 
% subset a-set = [x x X ...]. 
% 


% ndof: Total number of DOF. Set is labeled "nset". 
% aset: Retained DOF (proper subset of [1:1:ndof]) 


nset = [1 :ndof]; 
for icnt = 1 : length(aset); 

indices(icnt) = find(nset == aset(icnt)), 
end 


bool = ones(size(nset)); 


bool(indices) = zeros(size(indices)); 
oset = nset(find(bool>0)); 
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function[k,m]=fspringmass(conn,mspec,fix) 


% 
% 
% 
% 
% 
% 
% 
% 
% 
% 
% 
% 
% 
% 
% 
% 
% 
% 
% 
% 
% 
% 
% 
% 
% 
% 
% 


This function assembles the stiffness [k] and mass [m] matrices for an assemblage 
of springs. The connectivity is established in the array ‘conn’, the boundary 
conditions (fixity) are specified in the array 'fix'. The mass is specified in the 
vector mspec. 


For example, to build the following system: 


m 2m m 
|--////--* --//[]--*¥==1|/==* 
k k 2k 
the connectivity array is as follows: 
[12 
23 
34 
3 4] 
The connectivity array is therefore of dimension: 
(number of springs)*(2). 


The fixity array for the above system is 
fix=[0111] 


The fixity array is therefore of dimension: 
(1)*(number of nodes). 


The mass distribution is given by mspec: 
mspec = [1 2 1] 


The mspec array is of dimension: 
(1)*number of dynamic (free) DOF 


kspring=[1 -1;-1 1]; 
ksize=length(fix); 
k=zeros(ksize,ksize); 
m=zeros(length(mspec)); 


% assemble stiffness matrix: 


for spring = 1:size(conn, 1); 


index1=conn(spring,1); 
index2=conn(spring,2); 
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k(index 1 ,index 1)=k(index1,index1)+kspring(1,1); 
k(index 1 ,index2)=k(index 1 ,index2)+kspring(1,2); 
k(index?2, index 1 )=k(index2,index1)+kspring(2, 1); 
k(index2,index2)=k(index2,index2)+kspring(2,2); 
end: . 


% zero rows and columns corresponding to boundary conditions 
indices _to_keep = find(fix); 
if length(indices_to_keep) ~= 0; 
k =k(indices_to_keep,indices_to_keep); 
end 
%assemble mass matrix: 
m=zeros(size(k)); 
for i=1:length(m); 


m(i,i)=mspec(i); 
end 
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